

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%  NOTES ON COMPUTATION: THE DEFAULT BEHAVIOR OF THIS FILE IS TO COMPUTE ESTIMATES ONCE USING
%%%%%%%  THE KNOWN GLOBAL SOLUTION AS THE START POINT.  THIS BEHAVIOR CAN BE ALTERED (SAY FOR A
%%%%%%%  DIFFERENT DATASET OR FOR A DIFFERENT SPECIFICATION) BY CHANGING THE VALUE OF "AlreadyDone"
%%%%%%%  FROM "1" TO "0" INSTEAD, WHICH WILL TRIGGER A LONGER SOLUTION WITH MULTIPLE RESTARTS AND
%%%%%%%  VARIOUS OPTIMIZATION METHODS, AS CODED IN THE "else" PART OF THE LOOP BELOW.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


AlreadyDone = 1 ;
if AlreadyDone==1
    disp('Computing Plots with the last finalized point estimates...')
    load('PointEstimates_MSM7.mat')

    ThetaLEstimatorMSM(EPointEst,AchievementData,LPointEst.knotvec_c,[],LPointEst.params,1,0,0,[])
else
%%
    % poolobj    = parpool('local',96,'IdleTimeout',180) ;
    poolobj    = parpool('local',24,'IdleTimeout',180) ;
    % poolobj    = parpool('local') ;

%%


    timestamp1 = clock ;
    datetime(now,'ConvertFrom','datenum')  %#ok
    disp('***')
    disp('***')
    disp('***')
    disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
    disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
    disp('%%%%                            NOW COMPUTING STAGE-2 FIRST-PASS POINT ESTIMATES USING PATTERN SEARCH...                             %%%%')
    disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
    disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
    disp('***')
    disp('***')
    disp('***')
    tic
    LPointEst      = ThetaLEstimatorMSM(EPointEst,AchievementData,[],[],[],1,0)  %#ok
    ssrPreliminary = LPointEst.ssrvalc ;
    toc

    %%

    keepnum = 200 ;
    tic
    disp('***')
    disp('***')
    disp('***')
    disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
    disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
    disp('%%%%  NOW COMPUTING STAGE-2 MULTI-METHOD POINT ESTIMATES USING GENETIC ALGORITHM, ADAPTIVE START POINT SEARCH, AND PATTERN SEARCH... %%%%')
    disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
    disp('%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%')
    disp('***')
    disp('***')
    disp('***')
    cutoffs   = EPointEst.cutoffs ;
    knotvec_c = LPointEst.knotvec_c;
    pi_c      = LPointEst.params ;
    P         = ThetaLEstimatorMSM(EPointEst,AchievementData,knotvec_c,[],pi_c,0,0,1,[]) ;

    BSstate   = 10 ;  %%%%This Bootstrap-state switch value indicates that
    %%%%we are currently doing preliminary point estimates via genetic algorithm
    if length(knotvec_c)==11
        load('BootstrapStartValues.mat')
        pigrid = nan(length(Solutions{1,1}),length(Solutions)) ;
        for ii=1:length(Solutions)
            pigrid(:,ii) = Solutions{ii,1} ;
        end
        temp = size(pi_c) ;
        if temp(2) == 1
            pigrid = [pi_c pigrid] ;
        else
            pigrid = [pi_c' pigrid] ;
        end
    else
        pigrid = pi_c ;
    end
    [pigrid,ssrlog,survivorfrac] = ASPS(pigrid,knotvec_c,P,BSstate,keepnum,0) ;
    survivorfrac  %#ok

    GMMRound       = 1 ;
    stopflagPhase1 = 0 ;
    while stopflagPhase1==0 && GMMRound<=5
        disp('%%%%GAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGA%%%%%')
        disp(strcat('%%%% SIMULATED GMM ESTIMATION, ROUND_',num2str(GMMRound),'_VIA GENETIC ALGORITHM...'))
        disp('%%%%GAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGA%%%%%')
        LPointEsttemp = ThetaLEstimatorMSM(EPointEst,AchievementData,knotvec_c,[],[],0,BSstate,0,pigrid') ;

        if LPointEsttemp.exitflagc==4 || LPointEsttemp.exitflagc==5
            disp('%%%%%%%%%%%%%%%%%%%%PHASE_2_PHASE_2_PHASE_2_PHASE_2_PHASE_2_PHASE_2_PHASE_2_%%%%%%%%%%%%%%%%%%%%%')
            disp('%%%%%%%%%%%%%%%%%%%%     INITIATING TERMINAL PHASE 2 OF ESTIMATION...       %%%%%%%%%%%%%%%%%%%%%')
            disp('%%%%%%%%%%%%%%%%%%%%PHASE_2_PHASE_2_PHASE_2_PHASE_2_PHASE_2_PHASE_2_PHASE_2_%%%%%%%%%%%%%%%%%%%%%')
            [pigrid,ssrlog,survivorfrac] = ASPS([LPointEsttemp.params; LPointEsttemp.population]',knotvec_c,P,BSstate,keepnum,1) ;  %#ok
            survivorfrac  %#ok
            ssrRelativeRange = (max(ssrlog)-min(ssrlog))/min(ssrlog)   %#ok

            BSstate(1) = 11 ;  %%%%This Bootstrap-state switch value indicates that
            %%%%we are currently doing final point estimates via genetic algorithm
            disp('%%%%GAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGA%%%%%')
            disp('%%%%       SIMULATED GMM ESTIMATION, TERMINAL ROUND...    %%%%%')
            disp('%%%%GAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGA%%%%%')
            LPointEsttemp = ThetaLEstimatorMSM(EPointEst,AchievementData,knotvec_c,[],[],0,BSstate,0,pigrid') ;
            ssrGAterminal = LPointEsttemp.ssrvalc ;

            disp('%%%%PSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPS%%%%%')
            disp(strcat('%%%% TERMINAL ROUND 1 OF SIMULATED GMM ESTIMATION VIA PATTERN SEARCH...     %%%%%'))
            disp('%%%%PSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPS%%%%%')
            BSstate(1)       = 0 ;  %%%%This Bootstrap-state switch value indicates that
            %%%%we are currently doing final point estimates via pattern search
            LPointEsttemp    = ThetaLEstimatorMSM(EPointEst,AchievementData,knotvec_c,[],LPointEsttemp.params',0,BSstate,0,pigrid') ;
            ssrPSpenultimate = LPointEsttemp.ssrvalc ;

            [pigrid,ssrlog,survivorfrac] = ASPS(LPointEsttemp.params,knotvec_c,P,BSstate,keepnum,1) ;
            survivorfrac  %#ok
            ssrRelativeRange = (max(ssrlog)-min(ssrlog))/min(ssrlog)   %#ok

            disp('%%%%PSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPS%%%%%')
            disp(strcat('%%%% TERMINAL ROUND 2 OF SIMULATED GMM ESTIMATION VIA PATTERN SEARCH...     %%%%%'))
            disp('%%%%PSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPS%%%%%')
            [~,indx]   = mink(ssrlog,1) ;
            pi_c0      = pigrid(:,indx) ;
            BSstate(1) = 0 ;  %%%%This Bootstrap-state switch value indicates that
            %%%%we are currently doing final point estimates via pattern search
            LPointEsttemp = ThetaLEstimatorMSM(EPointEst,AchievementData,knotvec_c,[],pi_c0,1,BSstate,0,pigrid') ;
            break
        end

        [pigrid,ssrlog,survivorfrac] = ASPS([LPointEsttemp.params; LPointEsttemp.population]',knotvec_c,P,BSstate,keepnum,1) ;
        survivorfrac  %#ok
        ssrRelativeRange = (max(ssrlog)-min(ssrlog))/min(ssrlog)   %#ok

        if survivorfrac>=0.90 || ssrRelativeRange<0.05 || GMMRound==5
            stopflagPhase1 = 1 ;
            BSstate(1) = 11 ;
            disp('%%%%GAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGA%%%%%')
            disp('%%%%       SIMULATED GMM ESTIMATION, TERMINAL ROUND...    %%%%%')
            disp('%%%%GAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGAGA%%%%%')
            LPointEsttemp = ThetaLEstimatorMSM(EPointEst,AchievementData,knotvec_c,[],[],0,BSstate,0,pigrid') ;
            ssrGAterminal = LPointEsttemp.ssrvalc ;

            disp('%%%%PSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPS%%%%%')
            disp(strcat('%%%% TERMINAL ROUND 1 OF SIMULATED GMM ESTIMATION VIA PATTERN SEARCH...     %%%%%'))
            disp('%%%%PSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPS%%%%%')
            BSstate(1)       = 0 ;
            LPointEsttemp    = ThetaLEstimatorMSM(EPointEst,AchievementData,knotvec_c,[],LPointEsttemp.params',0,BSstate,0,pigrid') ;
            ssrPSpenultimate = LPointEsttemp.ssrvalc ;

            [pigrid,ssrlog,survivorfrac] = ASPS(LPointEsttemp.params,knotvec_c,P,BSstate,keepnum,1) ;
            survivorfrac  %#ok
            ssrRelativeRange = (max(ssrlog)-min(ssrlog))/min(ssrlog)   %#ok

            disp('%%%%PSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPS%%%%%')
            disp(strcat('%%%% TERMINAL ROUND 2 OF SIMULATED GMM ESTIMATION VIA PATTERN SEARCH...     %%%%%'))
            disp('%%%%PSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPSPS%%%%%')
            [~,indx]      = mink(ssrlog,1) ;
            pi_c0         = pigrid(:,indx) ;
            BSstate(1)    = 0 ;
            LPointEsttemp = ThetaLEstimatorMSM(EPointEst,AchievementData,knotvec_c,[],pi_c0,1,BSstate,0,pigrid') ;
        end
        GMMRound = GMMRound + 1 ;
    end


    LPointEst                      = LPointEsttemp ;
    LPointEst.ssrvalcGAterminal    = ssrGAterminal ;
    LPointEst.ssrvalcPreliminary   = ssrPreliminary ;
    LPointEst.ssrvalcPSpenultimate = ssrPSpenultimate ;
    timestamp2                     = clock ;
    ElapsedMinutes                 = etime(timestamp2,timestamp1)/60 %#ok
end







